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Two-dimensional  model  of  the  intrinsic  point  defects  behaviour 
during  Cz  silicon  crystals  growth 

Anatolii  I.  Prostomolotov,  Natalia  A.  Verezub 
Institute  for  Problems  in  Mechanics  RAS,  117526  Moscow,  pr.  Vemadskogo  101,  Russia 


ABSTRACT 

Two-dimensional  mathematical  model  of  the  intrinsic  point  defects  recombination  during  Cz  growth  of  dislocation-free 
silicon  single  crystals  is  developed.  The  results  of  its  verification  are  compared  with  the  data  of  the  one-dimensional  model 
supposing  the  "fast"  vacancies  and  interstitials  recombination  near  the  liquid-solid  interface.  For  various  growth  conditions 
and  with  using  of  the  calculated  two-dimensional  temperature  fields  in  Cz  silicon  crystals  the  resulting  distributions  of  these 
intrinsic  point  defects  in  a crystal  are  analyzed. 

Keywords:  two-dimensional  (2D)  model;  Czochralski  (Cz)  crystal  growth,  calculation;  silicon;  intrinsic  point  defects 
(IPD);  vacancy;  interstitial;  recombination. 

1.  INTRODUCTION 

Diffusive  and  convective  transport  of  the  intrinsic  point  defects  (IPD)  and  its  recombination  in  a growing  crystal  near  a 
liquid-solid  interface  (LSI)  causes  the  main  influence  on  peculiarity  of  a defect  formation  in  dislocation-free  silicon  single 
crystals  (formation  of  voids  and  oxygen  particles,  congestions  of  interstitials)1,2.  At  present  the  model  of  the  defect 
formation  in  dislocation-free  single  crystals  working  effectively  is  based  on  the  one-dimensional  approximation  supposing 
the  "fast"  IPD  recombination  and  using  transformation  of  two  basic  transfer  equations  to  one  equation  by  replacement  of 
variables2. 

However  in  the  two-dimensional  case  it  is  more  reasonable  to  solve  the  initial  system  of  two  equations,  applying  the 
standard  Newton  algorithm  linearizing  two  equations  with  the  subsequent  approximation  and  the  numerical  solution  by  the 
finite  elements  method3.  In  addition  there  is  a possibility  of  calculation  in  regions  with  curved  boundaries  by  taking  into 
account  the  curved  LSI  shape  and  by  using  the  temperature  and  temperature  gradients  distributions  calculated  by  using  the 
global  thermal  Cz  model4.  In  given  work  this  approach  was  developed. 

For  verification  of  the  two-dimensional  algorithm  the  case  of  the  one-dimensional  (given  analytically)  temperature  field  in  a 
crystal  was  applied,  because  in  this  case  the  concentration  of  IPD  fields  are  known  on  basis  of  the  one-dimensional 
modeling2,5.  The  input  parameters  influence  on  the  residual  vacancies  and  interstitials  concentrations  and  on  the 
recombination  threshold  are  investigated.  In  particular  ^ value  is  studied  according  to  the  known  transition  criterion  from 
the  "interstitials"  growth  mode  to  the  "vacancy"  one  (here  £ = Vp/G  and  Vp  is  the  pulling  rate,  G is  the  axial  temperature 
gradient  at  LSI).  The  further  calculations  are  carried  out  for  a two-dimensional  temperature  field,  appropriate  to  growing  of 
silicon  single  crystals  of  150  mm  in  diameter.  Two-dimensional  fields  of  the  IPD  concentrations  are  calculated  at  various 
growth  stages  (30,  50  and  80%  of  length  of  the  cylindrical  grown  crystal  part)  and  for  various  pull  rates. 

2.  FORMULATION  OF  THE  PROBLEM 

The  crystals  growth  from  a melt  by  Cz  method  is  carried  out  in  the  high-temperature  hot  zone,  which  thermal  efficiency 
depends  on  a design  of  thermal  shields  and  properties  of  chosen  thermoinsulating  materials.  Besides,  hot  zone  should  satisfy 
technological  conditions  of  perfect  crystal  growth.  In  ref  [4]  the  conjugate  mathematical  model  for  the  investigation  of  heat 
transfer  processes  in  hot  zone  for  the  Cz  method  was  developed.  In  the  radiative-conductive  approximation  temperature 
fields  and  the  distribution  of  G for  consecutive  growth  stages  of  silicon  single  crystals  of  150  mm  in  diameter  were 
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investigated.  Calculations  of  heat  transfer  processes  were  carried  out  for  10  consecutive  growth  stages  depending  on 
volume  of  a cylindrical  part  of  a growing  crystal  (in  percentage):  the  beginning  of  a cylindrical  part,  then  the  stages 
corresponding  to  the  increase  of  crystal  volume  on  every  10%.  The  regularities  of  axial  temperature  gradients  and  £,  values 
change  at  the  LSI  were  obtained  which  can  be  used  for  theoretical  estimations  of  the  IPD  recombination  processes  in 
silicon.  In  fig.  1 the  frame  of  a hot  zone  and  isotherms  are  shown  for  the  50%-s  crystal  growing  stage. 

During  Cz  silicon  crystal  growth  the  vacancies  and  interstitials  transfer  is  realized  by  the  convective  diffusion.  We  neglect 


z.  cm 

130 


120  h 


10  20 


r,  cm 


Of 

s 

Q 

ta 

m 

s 

a 

b 

B 

B 

a 

B 

B 

B 

B 

TA 

0 

TA 

B 

B 

1 TA 

a 

Bl 

B 

B 

B 

B 

B 

B 

B 

B 

01 

a 

B 

ta 

B 

a 

TA 

B 

B 

B 

B 

B 

B 

B 

B 

a 

01 

Q 

ta 

B 

B 

a 

(Bl 

B 

B 

B 

B 

B 

B 

B 

B 

0 

K 

Ta 

S 

S 

a 

B 

B 

B 

B 

B 

B 

B 

B 

01 

1 0 

ta 

Q 

ta 

b 

B 

B 

B 

B 

B 

B 

B 

0 

TA 

B 

a 

S 

B 

B 

B 

B 

B 

B 

B 

0 

ta 

ta 

ta 

B 

a 

B 

B 

B 

B 

B 

B 

B 

0 

0 

ta 

g 

B 

b 

B 

a 

B 

B 

B 

B 

B 

B 

B 

Q 

B 

Bl 

a 

B 

B 

B 

B 

B 

B 

9 

9 

ft 

Q 

y 

y 

8 

6 

B 

B 

B 

s 

9 

ft 

ft 

ft 

ft! 

0 

0 

0 

0 

0 

ft 

0 

ft 

ft 

ft 

ft 

0 

0 

0 

0 

0 

0 

0 

g 

ft 

ft 

ft 

0 

0 

0 

0 

0 

0 

g 

ft 

ft 

01 

01 

0 

s 

0 

0 

ft 

0 

0 

* 

ss 

§ 

01 

9 

0 

0 

a 

0 

ft 

0 

0 

0 

0 

ft 

SI 

0 

0 

ft 

0 

0 

0 

a 

g 

g 

g 

ss 

89 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

S 

a 

1 

i 

1 

* 

0 

0 

g 

S 

0 

0 

| 

0 

| 

* 

1 

I 

i 

■s 

g; 

g 

i 

§ 

■■ 

B1 


Hn 


B, 


/B2 


B, 


Fig.2.  Scheme  of  the  calculation  region  (on  the 
right)  and  according  mesh  of  6-nodes  triangular 
finite  elements  (on  the  left). 

Here:  - symmetry  axis,  B2  - LSI,  B3  is  crystal 

surface,  B4  is  top  boundary  of  a recombination. 


Fig.  1.  Temperature  distribution  in  the  hot  zone  for  the 
50%-s  crystal  growing  stage:  1 - single  crystal,  2 - melt, 
3 - liquid-solid  interface,  4 - crucible,  5 - susceptor,  6 - 
resistive  heater,  7 - thermal  shields,  8 - water-cooled 
tank. 


the  thermodiffusion.  The  stationary  equations  system  for  calculation  of  the  vacancies  Cv(r,z)  and  interstitials  C{(r,z) 
concentrations  has  the  following  form: 

Vp^  = div(DvVCv)-co  (1), 

F dz 

Vp-^  = div(DiVCi)-C0  (2) 

F dz 

Coefficients  of  diffusions  (Dv,  Dj)  and  equilibrium  concentrations  of  vacancies  and  interstitials  (Cve,  Cie)  are  set  by  using  its 
values  (C™  at  the  melting  temperature  ( T ^ as  follows: 

Dv(T)  = Dvmexp(-EvD/kT+EvD/kTm),  Di(T)  = Dimexp(-EiD/kT+EiD  /kTm), 

Cve(T)  = CV7n  exp(-Ev  /kT  + Ev  /kTm),  C,e(T)=C,m  expf-E,  /kT  + E,  /kTm) 

Here  to  = Kvi(CvCi-CveCje),  the  recombination  factor  = Av](Dv  + D,)exp(-Erec/kT)  and  4,  = Anrcap  ( rcap - 
radius  and  Erec-  recombination  energy);  EvD,  EiD  - energies  of  vacancies  and  interstitials  diffusion  activation,  and  Ev,  E,  - 
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energies  of  vacancies  and  interstitials  equilibrium  activation.  In  these  formulas  the  constants  are  used:  silicon  melting 
temperature  Tm  =1683^T  and  Boltzman  constant  k = 8.625 xl(T5  eV/K. 

The  experimental  value  of  a transition  - %crit  =0.132  mm2/Kmin  from  vacancies  to  interstitials  growth  mode  corresponds  to 
the  following  parameters5: 

= 4.0xl0’5cm2/s,EvD  ^035  eV,  Dim  = 5.25x1 0~4 cm2/s,  EiD  - 0.20  eV;  (3) 

= 7.2 x 1014  cm'3,  Ev  = 4.5  eV,  Clm  = 5.2 x 1014 cm3,  Ei  = 4.6  eV. 

The  value  of  recombination  energy  barrier  Erec  is  not  known.  At  its  choice  the  satisfaction  to  “fast”  recombination  condition 
is  required.  This  condition  is  satisfied,  if  Erec  = 1.5  eV  is  set.  Radius  of  recombination  r^is  equal  to  3.0  x 10"8  cm. 

The  scheme  of  calculation  region  is  given  on  Fig.  2.  Accordingly  boundary  conditions  are  set  as  follows: 

SC  ^ 

= 0 at  Bi,  B4  boundaries,  C(/)  = at  B2,  B3  boundaries. 

Sn 

Here  C0)  corresponds  to  the  actual  values  ofCv , C,  and  are  appropriate  equilibrium  concentrations.  For  calculations 
a geometry  and  sizes  of  a recombination  region  were  chosen  as  following  values:  H = 1 15.2  cm  and  Hr  = Ho  = 103.7  cm, 
according  to  the  size  of  global  thermal  Cz  model. 

The  stationary  equations  (1)  - (2)  are  solved  iteratavely  by  a finite  element  method  (FEM)  in  the  Galerkin’s  formulation. 
The  approximation  of  unknown  concentrations  Cv  and  Q is  represented  as  decomposition  on  square-law  basic  functions  on 
6-nodes  triangular  elements  (fig.  2).  Preliminary  a residual  for  each  equation  is  linearizated  by  Newton’s  method.  FEM 
equations  system  is  solved  by  a frontal  method.  As  a whole  the  algorithm  of  the  solution  of  the  discrete  FEM  equations, 
earlier  developed  for  global  thermal  model  is  kept4. 

3.  COMPARISON  WITH  ONE-DIMENSIONAL  MODEL 

For  the  test  based  on  the  data  of  one-dimensional  recombination  model  the  temperature  field  in  a crystal  was  set  as  only 
depending  upon  z: 

1 _ 1 Gz 

T(z)  ~Tm+  Tl  ’ 

SC® 

and  at  a lateral  crystal  surface  the  boundary  condition  was  set  as  follows: = 0 . 

Sn 

Here  G is  axial  gradient  of  temperature  at  LSI.  In  this  case  calculation  results  on  based  of  2D-model  should  correspond  to 
one-dimensional  data2,5. 

The  comparison  of  Vp/G  critical  value  at  G changes  in  range  of  30-300  K/cm  and  at  Vp  pull  rate  variation  near  (v-i)  - 
transition  threshold  shows,  that  the  transition  occurs  at  constant  value  of  (V  p/G)crit  = 0.132  mm2/Kmin,  that  precisely 
corresponds  to  the  data5. 

In  a case  of  2D  model  a value  of  at  Tm  was  equal  « 2x10" 19  cm3/s  and  for  this  value  a “fast”  recombination  condition 
of  CveCie  = CvCi  is  satisfied  with  rather  high  accuracy.  Also  a value  of  a residual  concentration  for  a vacancies  mode 

(G  = 120  K/cm,  Vp  = 1.61  mm/min)  was  equal  to  0*  = 3.9xl012  cm'3.  According  to  ref  [5]  this  value  is  a little  lower:  = 

2.8x1 012  cm'3.  The  small  divergence  can  be  caused  by  distinction  in  approximation  of  the  basic  equations,  because  in  ref 
[2,5]  the  "fast"  recombination  condition  is  used  for  transformation  of  (1)  - (2)  system  to  one  equation.  Opposite,  in  given 
2D-  model  the  initial  (l)-(2)  system  is  solved  and  the  "fast”  recombination  condition  follows  from  the  calculations. 

Thus,  it  is  possible  to  consider,  that  for  an  one-dimensional  temperature  field  the  received  results  rather  well  correspond  to 
lD-data. 


4.  IPD  BEHAVIOUR  IN  TWO-DIMENSIONAL  THERMAL  FIELD 

In  a case  of  2D  thermal  field  the  boundary  Cv,  C,  concentrations  are  set  as  Cv  = Cve  and  Q = C)e  at  a lateral  crystal  surface. 
Input  2D  thermal  fields  are  set  according  to  results  of  calculations  by  a global  thermal  model4  for  different  crystal  growth 
stages.  In  this  case  the  radial  inhomogeneities  of  an  axial  temperature  gradient  G at  LSI  plays  the  important  role.  Therefore 
determining  its  are  two  values  of  ^ = Vp/G0  - at  an  axis  and  = Vp/GR  - at  a lateral  crystal  surface.  Its  ratio  with 
determines  the  recombination  results  and  a growth  mode:  vacancy  (v),  interstitials  (i)  or  mixed  (v-i).  For  the  latter  case 
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a so-called  OSF-ring  is  created  in  a crystal.  Radius  of  this  ring  is  determined  experimentally  and  this  characteristic  plays  an 
important  role  at  IPD  distribution  verification. 


Fig.3.  Isolines  of  the  IPD  concentration  difference  (C  = Q - Cv  [cm°])  for  various  crystal  growth  stages  at  the  interstitial 
growth  mode:  a - 30%  body,  Vp  = 0.39  mm/min;  b - 50%,  Vp  = 0.37  mm/min. 

According  to  (3)  the  parameters  of  2D  model  correspond  to  4^  =0.132  mm2/Kmin.  However  for  2D  thermal  field  a value 

of  4cm  requires  an  additional  check.  For  this  purpose  the  appropriate  parametrical  calculations  were  carried  out,  in  which  the 
pull  rate  Vp  and  values  of  temperature  gradients  (G0  and  Gr)  were  varied  at  various  crystal  growth  stages.  In  a number  of 
calculations  (50-  and  80%-s'  growth  stages)  the  value  of  4cm  was  checked  by  means  of  Vp  change  for  conditions  according 
to  an  identical  temperature  field  in  a crystal,  i.e.  without  the  appropriate  recalculation  of  a thermal  field  on  global  thermal 
model. 

For  30%  and  50%  growth  stages  the  appropriate  values  Vp  = 0.39  and  0.37  mm/min  and  the  axial  temperature  gradients  are 
35  and  30  K/cm  at  an  axis,  51  and  42  K/cm  at  a crystal  edge.  At  30  % stage:  4o = 0.1 1 1,  4r  = 0.077,  and  at  50  %:  4o  = 0.126, 
4r  = 0.088  mm2/Kmin,  i.e.  at  B4  boundary  is  realized  the  interstitials  mode  (fig.  3a,b),  because  4 < &*• 


Fig.4.  Isolines  of  the  IPD  concentration  difference  for  the  50%-s  growing  stage  and  pull  rates  according  to  the 
mixed  (v-i)  and  vacancy  (v)  growth  modes:  a - Vp  = 0.42  mm/min;  b - 0.47  mm/min. 
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All  calculations  were  carried  out  at  the  initial  conditions  according  to  equilibrium  concentration  values  of  vacancies  and 
interstitials  in  calculation  region.  For  the  chosen  parameters  (3)  the  equilibrium  vacancies  concentration  exceed  the 
equilibrium  interstitials  concentration.  Therefore  for  initial  iterations  the  prevalence  of  vacancies  concentration  is  observed. 
For  the  subsequent  iterations  an  action  of  the  interstitials  source  in  a region  of  the  greatest  temperature  gradients  (near  LSI 
and  a crystal  edge)  leads  to  a visible  effect.  It  results  in  characteristic  elliptical  shapes  of  concentration  isolines  C = Q-  Cv. 
On  fig.  4a,b  the  vacancies  and  interstitials  distributions  are  shown  at  pull  rates,  which  are  a little  larger,  than  in  the  previous 
case  for  the  50%-s  growth  stage:  Vp  = 0.42  and  0.47  mm/min.  In  this  case  thermal  data  corresponds  to  former  value  Vp  - 
0.37  mm/min,  and  the  change  of  pull  rate  is  taken  into  account  only  in  the  recombination  equations  (l)-(2). 


Fig.  5.  Isolines  of  the  IPD  concentration 
difference  for  various  crystal  growth  stages  at 
Vp  = 0.50  mm/min  corresponding  to  the 
following  growth  modes: 
a - transition  from  the  mixed  (v-z)  to  the 
vacancy  (v)  mode  at  the  50%-s  stage; 
b - formation  of  the  mixed  (v-z)  mode  at  the 
80%-s  stage. 


The  growth  of  Vp  increases  = 0.143,  = 0.100  and  £*>  = 0.159,  = 0.1 12  mm2/Kmin,  accordingly.  In  particular  at  Vp 

= 0.42  mm/min  is  created  a mixed  v-i  mode,  which  turns  into  a vacancies  mode  at  increasing  of  pull  rate  till  0.47  mm/min. 
Let’s  notice,  that  remains  smaller  than  2^*. 

Strictly  speaking,  a pull  rate  change  causes  the  changes  of  a temperature  field  in  a crystal.  These  changes  my  be  accounted 
for  50%  growth  stage  at  Vp  = 0.50  mm/min.  The  axial  temperature  gradient  increases  then  to  32  K/cm  et  the  center  and  to 
46  K/cm  at  the  crystal  edge.  It  gives  next  values:  £*>  = 0.159  and  = 0.110  mm2/Kmin.  The  interstitials  region  extends 
essentially  such  that  the  vacancies  distribution  is  divided  from  a top  boundary  (fig.  5a). 
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Similar  changes  occur  at  increase  of  pull  rate  for  the  80%-s  growth  stage.  At  Vp  = 0.35  mm/min  the  appropriate  axial 
gradients  at  the  axis  and  at  the  crystal  edge  are  equal  G0  = 22  and  GR  = 48  K/cm,  and  ^ = 0. 1 10  and  = 0.073  nrni2/Kmin. 
It  may  be  a confirmation  of  an  interstitials  growth  mode. 

However  at  larger  pull  rate  - Vp  = 0.50  mm/min,  the  values  % grow  so  that  there  is  mixed  v-i  mode  = 0.157,  = 0.104 

mm2/Kmin).  Respective  calculation  results  are  shown  in  fig.  5b. 

The  main  results  of  solution  of  this  stationary  recombination  problem  are  the  radial  distribution  of  residual  concentrations  of 
vacancies  Cv  (r)  and  interstitials  atoms  Q (r)  at  the  top  boundary  of  calculation  region  B4  (Fig.  2).  These  data  are  necessary 
for  further  plotting  of  summarizing  2D  distributions  of  residual  IPD  defects  in  the  crystal  (which  are  plotted  on  based  of 
calculated  data  for  10  and  more  growth  stages). 

Here  for  comparison  the  radial  output  distributions  Cv  (r)  and  Q (r)  on  the  top  boundary  of  recombination  region  are  given 
in  fig.  6 for  the  50%  growth  stage  and  two  pull  rate  values,  according  to  interstitials  (Vp  - 0.37  mm/min)  and  mixed  v-i  (Vp 
= 0.50  mm/min)  modes.  The  analysis  of  the  graphs  at  r = 0 shows,  that  in  the  interstitial  mode  Q exceeds  Cv  concentration 
essensially.  The  similar  graphs  for  mixed  (v-i)  mode  show  significant  growth  of  vacancies  concentration  and  practically 
complete  disappearance  of  interstitial  atoms  at  the  axis. 


cv,c, 


Fig.6.  Radial  distributions  of  Cv,  and  Q 
concentrations  on  the  top  boundary  (z  ==  H) 
corresponding  to  pull  rate  variations  at  the  50%-s 
crystal  growing  stage.  Here:  the  interstitial  (i) 
mode  is  at  Vp  — 0.37  mm/min,  and  the  mixed  (v-i) 
mode  arises  at  Vp  = 0.50  mm/min. 


5.  CONCLUSIONS 

Considered  2D  model  develops  the  previous  ID  model  of  IPD  behaviour  near  the  LSI  during  Cz  silicon  crystal  growth.  The 
assumption  of  the  "fast"  vacancies  and  interstitials  recombination  taken  into  account  at  a choice  of  parameters  by  using  of 
ID  data  was  confirmed  numerically.  The  values  of  the  residual  IPD  concentrations  (after  a recombination)  were  verified 
too. 

Using  2D  global  thermal  Cz  model  for  calculations  of  2D  thermal  growth  history  gives  the  possibility  of  IPD  behaviour 
modeling  adequately  at  various  changes  of  thermal  shields  in  a hot  zone  and  for  different  pull  rates.  In  particular  this  model 
gives  possibility  for  investigation  of  important  2D  effects  such  as  OSF  ring  in  dislocation-free  silicon  single  crystals  caused 
by  radial  temperature  gradients. 
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